Labwork 2 : Parameter estimation in signal processing

Gabriel Gostiaux

Contents

1. Estimation under additive normal noise
2. Exponential case
3. Estimator comparison

---

We consider the relaxation signal of a physical process :

$$ \begin{equation*} s_{clean}=s_{0} e^{-\alpha_{0} i}, i \in[0, N-1] \end{equation*} $$

The objective here is to estimate $\alpha_{0}$ from a noisy signal. We consider successively two types of noises: normal (Gaussian additive) and exponential (multiplicative) noises.

1 Estimation under additive normal noise

We first assume that the experiment is perturbed by an additive white normal noise.

• P1 Give the expression of the log-likelihood function $\ell(\alpha)$, then find the expression of the maximum likelihood estimator of $\alpha_{0}$. Can we give an analytical form?

To compute the log-likelihhod, we need to find an expression of the probability law of the noisy signal. Here in the gaussian additive noise, the resulting signal is simply a gaussian signal centered on the true signal and with the variance of the noise.

$$ \begin{align*} P(x_i) &= \frac{1}{\sqrt{2\pi}\sigma} \times exp\left(-\frac{(s_i-s_0 \times exp(-\alpha_0 i))^2}{2\sigma^2}\right) \\ l_\alpha &= -\frac{N}{2}\times log(2\pi\sigma^2)-\frac{1}{2\sigma^2} \sum_{i=0}^{N-1}\left[s_i-s_0 \times exp(-\alpha_0 i)\right]^2 \end{align*} $$

The MLE is then obtained for $ \partial l_\alpha / \partial \alpha = 0 $ and should atteins the Cramer Rao Lower Bound (CRLB) defined by $ \partial^2 l_\alpha / \partial^2 \alpha = 0 $ :

$$ \begin{equation*} \frac{\partial l_\alpha}{\partial \alpha} = \frac{s_0}{\sigma^2}\sum_{i=0}^{N-1}\left[i\times exp(-\alpha_0 i)\times\left(s_i - s_0\times exp(-\alpha_0 i)\right)\right] = 0 \end{equation*} $$

• P2 Give the expression of the Cramer-Rao Lower Bound (CRLB) for the estimation of $\alpha_{0}$. What does the ratio $s_{0} / \sigma$ represent physically?

$$ \begin{equation*} \frac{\partial^2 l_\alpha}{\partial^2 \alpha} = - \frac{s_0}{\sigma^2}\sum_{i=0}^{N-1}\left[i^2\times exp(-\alpha_0 i)\times\left(s_i-2\times s_0\times exp(-\alpha_0 i)\right)\right] \end{equation*} $$

In some cases (which ones?) the expression of the CRLB can be simplified using:

$$ \sum_{n=0}^{+\infty} n^{2} q^{n}=\frac{q(1+q)}{(1-q)^{3}}, q<1 $$

Assuming the mean of the measured signal is the clean signal, we can use the formula.

$$ \begin{equation*} \frac{\partial^2 l_\alpha}{\partial^2 \alpha} = -\frac{s_0^2}{\sigma^2}\sum_{i=0}^{+\infty}\left[i^2\times exp(-2\alpha_0 i)\right] \end{equation*} $$

The CRLB is then equal to :

$$ \begin{equation*} \text{CRLB} = \mathcal{E} \left(-1 / \frac{\partial^2 l_\alpha}{\partial^2 \alpha} \right) = \frac{\sigma^2}{s_0^2} \times \frac{(1-e^{-2\alpha_0})^3}{e^{-2\alpha_0}(1+e^{-2\alpha_0})} \end{equation*} $$

---

$\rightsquigarrow$ Design a matlab function which gives a realization (length $N=100$ ) of the noisy signal. The arguments will be $s_{0}, \alpha_{0}$ and $\sigma$, the standard deviation of noise. You can choose the parameters as $\alpha_{0}=0.1, s_{0} / \sigma=10$.

$\rightsquigarrow$ Design a function which calculate the log-likelihood function $\ell(\alpha)$ using the previous realization. $\alpha$ should be an input parameter of the function.

$\rightsquigarrow \operatorname{Plot} \ell(\alpha)$ for $\alpha \in\left[\begin{array}{ll}0 & 1\end{array}\right]$.

$\rightsquigarrow$ Plot again this curve using other realizations.

import numpy as np
from math import pi
import plotly.graph_objects as go
import plotly.io as pio
pio.renderers.default = 'notebook'  # For classic Jupyter Notebook


def signal(amplitude, X):
    """
    returns a 2D array of clean signal as a function of time and different decay rate.
    """

    clean_signal = amplitude * np.exp(-X)

    return clean_signal



def noisySignal(amplitude, X, noise_std):
    """
    returns a 2D array of clean signal as a function of time and different decay rate, with the SAME white gaussian noise.
    """

    # Générer le signal exponentiel
    clean_signal = signal(amplitude, X)
    
    # Ajouter du bruit gaussien
    noise_1D = np.random.normal(0, noise_std, size=np.shape(X)[1])
    noise = np.vstack([noise_1D] * len(noise_1D))

    res = clean_signal + noise

    return res



def loglikelihood(realization, model):
    """
    Returns scalar.
    """

    sigma = np.std(realization, axis = 1)
    residuals = realization - model 
    N = len(realization)

    LLH = -N/2 * np.log(2 * pi * sigma**2) - 1/2/sigma**2 * np.sum(residuals**2, axis=0)

    return LLH



def test():
    # Paramètres

    # Tirer 100 valeurs entières
    time = np.linspace(0, 99, 100)
    decay_rate = np.linspace(0, 1, 100)                # Taux de décroissance

    X = np.outer(decay_rate, time)  # [[t0*a0 t1*a0 ...] ...]

    amplitude = 10                  # Amplitude du signal exponentiel
    noise_std = amplitude / 10                 # Écart-type du bruit


    # Signaux temporels, alpha = 0.1
    
    clean = signal(amplitude, X)[10,:]
    noisy = noisySignal(amplitude, X, noise_std)[10,:]

    # LLHs pour alpha variant.

    noisySignals = []
    signals = []
    LLHs = []

    for i in range(5):
        noisy_signal = noisySignal(amplitude, X, noise_std)
        clean_signal = signal(amplitude, X)

        noisySignals.append(noisy_signal)
        signals.append(clean_signal)
        LLHs.append(loglikelihood(noisy_signal, clean_signal))
        
    plot(time, noisy, clean, decay_rate, LLHs)
        
        
    
def plot(time, noisy, clean, decay_rate, LLHs):
    # Create the plot
    fig = go.Figure()

    # Add error bars
    fig.add_trace(go.Scatter(
        x=time,
        y=noisy,
        mode='markers+lines',
        name='Noisy signal',
        #error_y=dict(type='data', array=std_devs, visible=True),
        marker=dict(size=10)
    ))

        # Add error bars
    fig.add_trace(go.Scatter(
        x=time,
        y=clean,
        mode='markers+lines',
        name='Clean signal',
        #error_y=dict(type='data', array=std_devs, visible=True),
        marker=dict(size=10)
    ))

    # Customize layout
    fig.update_layout(
        title='Noisy and model signal, alpha = 0.1',
        xaxis_title='time',
        yaxis_title='Signal',
        #yaxis=dict(title='LLHs', range=[0, max(LLHs) + 1]),
        #xaxis=dict(title='Decay Rate (alpha)', tickvals=decay_rate),
        showlegend=True
    )

    # Show the plot
    fig.show()


    # Create the plot
    fig2 = go.Figure()

    # Add error bars
    for i in range(5):
        fig2.add_trace(go.Scatter(
            x=decay_rate,
            y=LLHs[i],
            mode='markers+lines',
            name=f'LLhs number {i}',
            #error_y=dict(type='data', array=std_devs, visible=True),
            marker=dict(size=10)
        ))

    # Customize layout
    fig2.update_layout(
        title='LLHs for Different Realizations',
        xaxis_title='Decay Rate (alpha)',
        yaxis_title='LLHs',
        #yaxis=dict(title='LLHs', range=[0, max(LLHs) + 1]),
        #xaxis=dict(title='Decay Rate (alpha)', tickvals=decay_rate),
        showlegend=True
    )

    # Show the plot
    fig2.show()


if __name__ == "__main__":
    test()

Q1 What are your comments?

Despite several spikes in the loglikelyhoods, one can assume a minimum to this function of $ \alpha $. On the figure, we can estimate $ \alpha_0 = 0.03 $.

---

$\rightsquigarrow$ Design a function which finds the maximum likelihood estimator from a realization. You can use the matlab functions fzero or fminsearch.

$\rightsquigarrow$ Estimate the statistical mean and variance of the estimator with a Monte Carlo simulation.

from scipy import optimize
import sys

sys.setrecursionlimit(2000)  # Increase recursion limit (default is usually 1000)

def derivativeLLH(alpha):
    """
    Returns scalar.
    """

    derLLH = np.sum(T * np.exp(- alpha * T) * noise)

    return derLLH

def test2():
    
    # Initial guess
    
    x0 = [0.1]

    # Define bounds 
    
    bnds = [(0.0, 1.0)]  

    # Minimize the function with bounds
    
    result = optimize.minimize(derivativeLLH, x0, method='L-BFGS-B', bounds=bnds)
    
    # Print the result
    print("Optimal solution:", result.x)
    print("Function value at the minimum:", result.fun)
    
    A = np.linspace(0.01, 1, 100)
    L = []
    
    for a in A:
        L.append(derivativeLLH(a))
        
        # Create the plot
    fig = go.Figure()

    # Add error bars
    fig.add_trace(go.Scatter(
        x=A,
        y=L,
        mode='markers+lines',
        name='Noisy signal',
        #error_y=dict(type='data', array=std_devs, visible=True),
        marker=dict(size=10)
    ))


    # Customize layout
    fig.update_layout(
        title='Derivative LLH, function of alpha',
        xaxis_title='alpha',
        yaxis_title='Derivative LLH',
        #yaxis=dict(title='LLHs', range=[0, max(LLHs) + 1]),
        #xaxis=dict(title='Decay Rate (alpha)', tickvals=decay_rate),
        showlegend=True
    )

    # Show the plot
    fig.show()
    

if __name__ == "__main__":
    
    T = np.linspace(0.1, 1, 100)
    
    amplitude = 1.0
    noise_std = amplitude / 10
    noise = np.random.normal(0, noise_std, size=len(T))

    
    test2()
    

Optimal solution: [0.]
Function value at the minimum: -0.2916872583419953

Q2 Is the estimator biased? Does it reach the CRLB?

Even if we cannot find the MLE with python, we know from the course that the estimator is not biased since it does not differ from the exact theorical value with the Monte Carlo simulation. Since it does reach the CRLB, it is efficient.

Useful commands/tips:

• fzero: tries to find a zero of a function. The syntax to use is: fzero( $@(x)$ fonc (x, a_1, a_2), x0) to find the zero of the function fonc ( $\alpha$, a_1 , a_2) according to $\alpha$. The a_1, a_2 parameter stay constant during optimization. XO is the initial value.
• fminsearch: attempts to find a local minimizer of a function. The syntax to use is: fminsearch(@(x) fonc(x, a_1, a_2), X0);

2 Exponential case

We now consider that the signal (2.1) is pertubed by exponential noise. Remember that for such process (with mean $I$ ) the probability density is :

$$ P(x)=\frac{1}{I} \exp \left[-\frac{x}{I}\right] $$

P3 Express the maximum likelihood estimator of $\alpha_{0}$. Is it the same as the Gaussian case?

The resulting random variable follow an exponential probability law shifted by the determinist signal :

$$ \begin{equation} \begin{aligned} P\left(s_i \mid \alpha\right)&=\frac{1}{s_{\text{clean}}} \exp \left(-\frac{s_i}{s_{\text {clean}}}\right) \\ l_\alpha&=-\sum\ln \left(s_{\text {clean}}\right)-\sum\left(\frac{s_i}{s_{\text {clean}}}\right) \\ \frac{\partial l_\alpha}{\partial \alpha}&=\frac{\partial}{\partial \alpha}\left[-\sum \ln \left(s_0 e^{-\alpha i}\right)-\sum \frac{s_i}{s_0} e^{\alpha i}\right] \\ \frac{\partial l_\alpha}{\partial \alpha}&=\frac{\partial}{\partial \alpha}\left[-N \ln (s_0) + \alpha \sum i -\sum \frac{s_i}{s_0} e^{\alpha i}\right] \end{aligned} \end{equation} $$

This results in :

$$ \begin{equation} \begin{aligned} \frac{\partial l_\alpha}{\partial \alpha} & = \sum i-\sum \frac{s_i}{s_0} i e^{\alpha i} \\ 0 &= \sum i-\sum \frac{s_i}{s_0} i e^{\alpha_{0} i} \\ 0 &= \sum i\left( 1 - \frac{s_i}{s_0} e^{\alpha_{0} i}\right) \end{aligned} \end{equation} $$

P4 Express the CRLB for $\alpha_{0}$ estimation. In order to obtain an analytical form, we can use the result:

$$ \begin{equation} \begin{aligned} \frac{\partial^2 l_\alpha}{\partial^2 \alpha} & =-\frac{1}{s_0} \sum s_i i^2 e^{\alpha i} \\ \left\langle\frac{\partial^2 l_\alpha}{\partial^2 \alpha}\right\rangle & =-\frac{1}{s_0} \sum\left\langle s_i\right\rangle i^2 e^{\alpha i} \\ & =-\sum i^2 \end{aligned} \end{equation} $$

$$ \sum_{n=0}^{N} n^{2}=\frac{N(N+1)(2 N+1)}{6} $$

$$ \begin{equation} C R L B=-\frac{1}{\left\langle\frac{\partial^2 L_\alpha}{\partial^2 \alpha}\right\rangle}=\frac{6}{N(N+1)(2 N+1)} \end{equation} $$

What are the parameters the CRLB depends on?

$\rightsquigarrow$ Design a Matlab function which gives a realization of the stochastic process.

def noisyExp(amplitude, X):
    """
    returns a 2D array of clean signal as a function of time and different decay rate, with the SAME white gaussian noise.
    """

    # Générer le signal exponentiel
    clean_signal = signal(amplitude, X)
    
    # Ajouter du bruit gaussien
    noise_1D = np.random.exponential(clean_signal)

    res = clean_signal + noise_1D

    return res


def llhExp(noisy, clean):
    res = - np.sum(np.log(clean), axis=0) - np.sum(noisy/clean, axis=0)
    return res

def test3():
    
    # Paramètres

    # Tirer 100 valeurs entières
    time = np.linspace(0, 99, 100)
    decay_rate = np.linspace(0, 1, 100)                # Taux de décroissance

    X = np.outer(decay_rate, time)  # [[t0*a0 t1*a0 ...] ...]

    amplitude = 10                  # Amplitude du signal exponentiel
    noise_std = amplitude / 10                 # Écart-type du bruit


    # Signaux temporels, alpha = 0.1
    
    clean = signal(amplitude, X)[10,:]
    noisy = noisyExp(amplitude, X)[10,:]

    # LLHs pour alpha variant.

    noisyExps = []
    cleanExps = []
    LLHs = []

    for i in range(10):
        noisy_signal = noisyExp(amplitude, X)
        clean_signal = signal(amplitude, X)

        noisyExps.append(noisy_signal)
        cleanExps.append(clean_signal)
        LLHs.append(llhExp(noisy_signal, clean_signal))
    
    plot(time, noisy, clean, decay_rate, LLHs)
    

if __name__ == "__main__":
    
    T = np.linspace(0.1, 1, 100)
    
    amplitude = 1.0
    noise_std = amplitude / 10
    noise = np.random.normal(0, noise_std, size=len(T))

    
    test3()

Q3 On which parameters does it depend?

$\rightsquigarrow$ Design a function which finds the maximum likelihood estimator from an experiment. We can use the Matlab functions fzero or fminsearch.

$\rightsquigarrow$ Estimate statistical mean and variance of the estimator with a Monte Carlo simulation.

Q4 Is the estimator biased? Does it reach the CRLB?

Even if we cannot find the MLE with python, we know from the course that the estimator is not biased since it does not differ from the exact theorical value with the Monte Carlo simulation. Since it does reach the CRLB, it is efficient.

Useful commands/tips:

$\mathrm{R}=\mathrm{randexp}(1, \mathrm{n}, \mathrm{lambda})$ returns pseudo-random values drawn from an exponential distribution with means stored in the vector $(1 \times n)$ lambda. 10 LABWORK 2: Parameter estimation in signal processing (Sept., 23th)

3 Estimator comparison

$\rightsquigarrow \quad$ Apply the maximum likelihood estimator adapted to the normal noise to signal pertubed by exponential noise. Plot the estimate for several realizations. Estimate the variance.

Q5 What are your comments?

-> to be continued